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The  work  presented  in  this  paper  focuses  on  the  design  of  an  attitude  determination  and 
control  subsystem  (ADCS)  for  a  proximity  operation  and  imaging  satellite  mission.  The 
ARAPAIMA  (Application  for  Resident  Space  Object  Proximity  Analysis  and  IMAging) 
mission  is  carried  out  by  a  6  U  CubeSat  class  satellite  equipped  with  a  warm  gas  propulsion 
system.  The  propulsion  system  comprises  an  orbital  maneuvering  thruster  which  produces 
100  mN  and  a  set  of  16  reaction  control  system  (RCS)  thrusters,  of  25  mN  each,  installed 
in  pairs  that  generate  torques  about  each  of  the  satellite  body  axes.  The  thrust  of  the  RCS 
thrusters  can  be  modulated  over  the  entire  range  in  steps  of  1%  due  to  the  rapid  solenoid 
valve  actuation.  The  requirement  of  the  control  system  is  to  provide  pointing  control 
accuracy  of  1  arcmin  at  3  cr  in  the  desired  imaging  direction.  The  ADCS  employs  two  control 
laws.  One  control  law,  for  large  angle  maneuvers,  implements  eigenaxis  maneuvering  and 
the  other,  for  accurate  pointing,  is  implemented  with  PID  controllers  about  each  body 
axis. 

Simulations  performed  for  tracking  the  resident  space  object  flying  in  a  circular  orbit 
of  500  km  altitude  from  a  relative  orbit  of  250  m  are  used  to  test  the  performance  of 
the  ADCS.  A  comparison  between  a  “traditional”  ADCS  system  with  reaction  wheels 
and  RCS  thrusters  for  their  off-loading,  and  a  system  with  only  RCS  thrusters  has  been 
performed.  The  accuracy  of  the  pointing  is  comparable,  the  pointing  performance  as  well 
as  the  propellant  usage  for  both  options  are  further  discussed  in  the  paper.  The  integration 
of  both  controllers  within  the  ADCS  and  the  switching  are  also  discussed. 


Nomenclature 


Cm  Moment  coefficient 

Isp  Propellant  specific  impulse 

Jp  Solar  panel  moment  of  inertia  matrix 

Jrw  Flywheel  moment  of  inertia 

Km  Reaction  wheel  motor  torque  constant 

Kn  Knudsen  number 

Msi  Propellant  slosh  disturbance  moment 

N  Atomic  nitrogen 
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O  Atomic  oxygen 

Rm  Reaction  wheel  motor  resistance 

S  Maximum  projected  satellite  surface  area 
Ss  Sunlit  surface  area 

T  RCS  thrust 

a  Angle  of  attack 

/3  Sideslip  angle 

Be  Earth’s  magnetic  field  vector  at  satellite  location 

Fb  Force  vector  acting  on  the  body-frame 

Fd  Disturbance  force  vector  acting  on  the  body-frame 

Fcmd  Controller  Force  Command 

Hruj  Reaction  wheel  angular  momentum 

Jrw  Reaction  wheel  moment  of  inertia  matrix 

J  Moment  of  inertia  matrix 

Mb  Moment  vector  acting  on  the  body-frame 

Md  Disturbance  moment  vector  acting  on  the  body-frame 

Maero  Aerodynamic  disturbance  moment 

Mcmd  Controller  moment  command 

Mgg  Gravity  gradient  disturbance  moment 

Mp  Solar  panel  flexible  mode  disturbance  moment 

-M-rmm  Magnetic  disturbance  moment 

Msrp  Solar  radiation  pressure  disturbance  moment 

FCmd  Controller  thrust  command  at  time 

Vqo  Freestream  velocity  vector 

fir™  Reaction  wheel  angular  velocity  vector 

<I>  Solar  intensity  vector 

Q  Angular  velocity  imaginary  quaternion 

q  Unit  quaternion/Euler- Rodrigues  parameters 

Vb  Velocity  vector  of  the  body-frame 

e  Euler-axis  unit  vector 

r  Satellite  position  unit  vector  in  the  ECI  frame 

a?  Angular  velocity  vector  of  the  body-frame 

0P  Solar  panel  position  vector 

rrirmm  Residual  magnetic  dipole  or  moment 
qv  Vector  quaternion  component 

r  Satellite’s  position  vector  in  the  ECI  frame 

A  Atmospheric  particle  mean  free  path 

/d  Earth’s  gravitational  parameter 

l de  Earth’s  angular  velocity  vector 

ujrw  Reaction  wheel  angular  speed 

</>  Euler  angle  parameter 

p  Mass  density 

6 si  Mechanical  slosh  model  angular  displacement 

£  Reflectance  value  of  the  satellite 

b  Solar  panel  joint  damping  coefficient 

c  Speed  of  light  in  a  vacuum 

fs  Average  solar  constant 

frcs  RCS  thruster  operating  frequency 

ge  Gravity  at  Earth’s  surface 

hrw  Reaction  wheel  nominal  angular  momentum 

hsi:n  Mechanical  slosh  model  angular  momentum 

k  Solar  panel  joint  spring  constant 

ki:rw  Reaction  wheel  controller  integral  gain 

kp,rw  Reaction  wheel  controller  proportional  gain 

l  Reference  length  (shortest  satellite  body  dimension) 
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ls  Solar  radiation  pressure  moment  arm 
lrcs  RCS  thruster  moment  arm 

m  Satellite  rigid-body  mass 

m\)^rw  Propellant  mass  per  reaction  wheel  off-load  maneuvers 

niprop  Propellant  mass 

msi  Mechanical  slosh  model  mass 

n  Number  density 

qo  Scalar  quaternion  component 

s  Laplace  transform  complex  parameter 

xsi:Tl  Mechanical  slosh  model  linear  displacement 

I.  Introduction 

With  the  increasing  interest  on  low  budget  space  missions,  academia  has  found  a  niche  in  the  space 
industry.  Going  further  than  proposals  and  designs,  many  universities  have  flown  satellites  as  parts  of  military 
and  governmental  initiatives  across  the  world.  These  satellites,  better  known  as  CubeSats  (minuaturized 
satellites  measured  in  1  liter  volume  units  or  1  U)  ,  can  be  developed  using  commercial  off-the-shelf  (COTS) 
products  which  considerably  lower  the  overall  cost  of  missions.  Recently,  CubeSat  missions  have  grown 
more  ambitious  in  nature  helped  along  by  the  current  trend  on  miniaturized  high  performance  optics,  power 
systems,  and  other  electronics.  These  missions  now  blur  the  line  which  distinguishes  professional  projects 
from  those  previously  regarded  as  overreaching. 

While  encouraging,  this  upsurge  brings  a  set  of  issues  to  address.  As  popularity  and  feasibility  continue 
to  rise  so  does  the  number  of  satellites  launched  each  year.  Each  of  these  contribute  to  the  amount  of  debris 
in  low  Earth  orbit  (LEO)  (altitudes  ranging  from  160  km  to  2000  km).  Although  this  specific  issue  has 
recently  been  the  focus  of  some  publicity  no  final  solution  currently  exists. 

A.  Arapaima  Mission  Overview 

As  a  precursory  component  to  the  reduction  of  orbital  debris,  the  ARAPAIMA  mission  proposes  a  reconnais¬ 
sance  approach  to  perform  visible,  infrared  (IR),  and  3D  imaging  of  Resident  Space  Objects  (RSOs)  without 
a  priori  knowledge  of  their  shape  or  attitude.  This  process  follows  a  set  of  autonomous  approach  and  close 
proximity  maneuvers  carried  out  with  sufficient  accuracy  to  allow  rendezvous  and  docking  maneuvers  with 
the  RSO. 

The  mission  is  carried  out  by  a  12  kg  6  U  (36.6  x  23.9  x  27.97)  cm  CubeSat  (with  an  imaging  array 
consisting  of  an  IR  camera,  a  miniature  laser  rangefinder,  and  a  visible  light  monochrome  camera  arranged 
such  that  their  respective  imaging  directions  are  parallel  to  each  other  as  shown  in  figure  1. 

To  perform  orbital  approach  maneuvers,  the  satellite  is  equipped  with  a  difluoroethane  warm  gas  propul¬ 
sion  system  operated  via  rapid  solenoid  valve  actuation  and  miniaturized  2D  nozzle  sets  attached  to  the 
satellite  body.  The  system  is  comprised  of  16  RCS  thrusters  setup  in  pairs,  each  one  capable  of  producing 
up  to  25  mN  of  thrust.  The  thrusters  are  arranged  such  that  4  nozzles  are  parallel  to  each  other  in  the  x  and 
z  directions  on  the  body  frame  (figure  2).  These  4  nozzle  clusters  make  up  the  orbital  maneuvering  thrusters 
(OMTs)  and  provide  a  total  100  mN  of  thrust. 

Reaction  attitude  control  is  performed  by  operating  the  thrusters  in  pairs;  these  are  positioned  such  that 
a  pure  moment  can  be  exerted  about  each  body  axis.  The  number  of  thrusters  and  their  placement  provide 
a  full  layer  of  redundancy  as  the  RCS  system  is  able  to  operate  with  4  nozzle  pairs  (pairs  A  B  C  D  or  E  F 
G  H  in  figure  2). 

The  satellite’s  attitude  determination  and  control  system  (ADCS)  uses  a  combination  of  startracker 
(STR)  and  inertial  measurement  unit  (IMU)  sensor  readings  to  achieve  a  pointing  control  accuracy  of  1 
arcmin  at  3cr  in  the  desired  direction  of  the  imaging  array  (figure  1).  The  startracker  has  a  boresight 
accuracy  of  6  arcsec  and  a  40  arcsec  roll  axis  accuracy.  The  IMU  includes  angular  rate  gyro  (±150%  range) 
and  accelerometer  (±2  g  range)  triads  together  with  tri-axial  magnetometer  sensors  (±1.9  gauss).  As  a 
secondary  attitude  determination  system,  the  satellite  will  also  carry  photodiode  sun  sensors.  These  are 
used  during  operational  modes  in  which  the  startracker  is  unable  to  provide  an  attitude  solution  such  as 
detumbling  maneuvers.  An  onboard  GPS  module  is  utilized  to  track  the  position  of  the  satellite  during 
orbital  maneuvers  and  complement  the  communications  array. 
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MC:  Monochrome  camera  frame 
STR:  Star-tracker  frame 
©  Center  of  mass 


Figure  1:  Imaging  array  placement  on  the  satellite  body.  The  x-axis  on  the 
satellite’s  body-fixed  frame  is  defined  as  being  parallel  to  the  imaging  direction. 


®  Center  of  mass 
@  Solar  panel  array 


satellite  body.  Pairs  C,  D,  G,  and 
H  lie  on  the  +Yb  face. 


Figure  2:  RCS  thruster  placement  on  satellite  body-frame. 
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II.  Dynamic  Modeling 


The  satellite  dynamics  serving  as  the  process  or  plant  for  the  system  were  modeled  as  a  combination 
of  attitude  kinematics  and  orbital  mechanics  for  a  rigid  body  with  6  degrees-of- freedom  (DOF).  These 
are  further  augmented  by  a  set  of  environmental  models  and  conditions  which  contribute  to  the  external 
disturbances  acting  on  the  satellite  body.  Sources  of  internal  disturbances  are  also  considered. 

A.  Rigid-Body  Dynamics 

The  satellite  body  is  treated  as  a  rigid  body  with  constant  mass  and  moment  of  inertia.  The  forces  and 
moments  acting  on  the  satellite’s  body  fixed  frame  are  given  by  Eqs.  (1)  and  (2)  with  the  body-fixed  frame 
defined  as  shown  in  figure  2. 

Fb  =  m(Vb  +  u  x  Vb)  +  Fd  + 

Fcrnd  (i) 

Mb  =  J Co  x  (J uj)  +  Mf]  +  Mcrri(i  (2) 

Changes  throughout  the  mission  to  the  satellite’s  configuration  such  as  when  appendages  (solar  panels, 
antennas,  etc.)  are  deployed  are  modeled  as  impulsive  changes  in  its  moment  of  inertia  by  corresponding 
moment  disturbances  .  Additionally,  for  the  purpose  of  attitude  control  simulations  the  mass  flow  rate  of 
the  propellant  is  assumed  to  be  small  enough  for  the  mass  to  be  considered  constant  throughout  imaging 
maneuvers. 

The  attitude  kinematics  are  described  in  terms  of  Euler-Rodrigues  parameters  or  rotation  quaternions. 
All  quaternions  referenced  henceforth  are  considered  to  be  normalized  unit-quaternions  and  are  defined  as 
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where  e  is  an  arbitrary  unit  column  vector  and  <f>  is  an  arbitrary  angle  through  which  a  3-dimensional  frame 
is  rotated  about  e  as  given  by  Euler’s  formula1 


R(e,  =  cos ((j))I  +  (1  —  cos (<j)))eeT  +  cos  0  (4) 

The  quaternion  kinematic  equation  in  terms  of  the  inertial-referenced  body  angular  velocity  is  then  given 


by 
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Or  in  matrix  notation 
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Equations  (1)  and  (2)  along  with  Eq.  (5)  make  up  the  full  dynamics  of  the  model  and  allow  us  to  determine 
the  satellite’s  state  at  any  time  t  based  on  known  initial  conditions.  All  perturbations  and  inputs  to  the 
system  are  included  directly  into  Eqs.  (1)  and  (2)  through  the  Fd  and  Md  terms  respectively  and  are 
propagated  through  the  integration  of  u>. 
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The  error  quaternion  of  the  satellite  is  defined  as  the  quaternion  which  describes  the  rotation  from  the 
current  quaternion  state  to  the  desired  or  commanded  quaternion.  It  is  defined  by  quaternion  algebra  as 

qe  =  (qo,e )  =  (8) 

\Qv,eJ 

The  pointing  error  is  defined  as  the  angle  from  the  body  x-axis  (imaging  direction)  to  the  desired  x-axis 
as  given  by  the  error  quaternion.  It  can  be  easily  calculated  as 

V’e  =  (9o2,e  +  9b-de-de)  (9) 


B.  Environmental  Models 

To  account  for  environmental  conditions  and  scenarios  a  set  of  models  are  set  parallel  with  the  dynamic 
equations  discussed  above.  These  shape  the  moments  and  forces  acting  on  the  system  to  mirror  the  scenario 
of  a  satellite  in  LEO  at  approximately  500  km  from  the  Earth’s  surface.  At  this  altitude,  the  effects  of 
Earth’s  gravitational  and  magnetic  field  as  well  as  its  atmosphere  are  non-negligible.  As  such,  they  form  the 
majority  of  the  plant  or  process  under  consideration. 

At  this  point  in  time  it  would  also  be 
useful  to  specify  the  reference  frames  utilized 
throughout  the  different  models.  The  body 
frame  has  already  been  defined  as  seen  in  fig¬ 
ures  1  and  2.  The  inertial  frame  is  defined 
as  the  ECI  (Earth- Centered  Inertial)  as  de¬ 
scribed  by  the  J2000  system3  where  the  x-axis 
is  aligned  with  the  mean  vernal  equinox  and 
the  z-axis  is  aligned  with  the  uje  vector  shown 
in  figure  3.  This  frame  is  used  for  all  inter¬ 
nal  calculations  throughout  the  models  and 
it  serves  as  a  reference  point  for  transforma¬ 
tions  between  different  frames.  Another  use¬ 
ful  frame  is  the  ECEF  (Earth- Centered  Earth- 
Fixed)  where  the  x-axis  points  to  the  (0°,  0°) 
point  on  the  Earth’s  graticule  and  the  z-axis 
is  aligned  with  the  uje  vector.  It  should  be 
noted,  however,  that  because  ECEF  frame  ro¬ 
tates  with  the  Earth,  the  z-axis  does  not  lie 
perfectly  along  uje  due  to  polar  motion. 

Finally,  some  environmental  models  pro¬ 
vide  information  on  the  North-East-Down 
(NED)  frame,  a  type  of  noninertial,  flat  Earth 
(local  tangent  plane,  LTP)  reference  system 
usually  used  for  aircraft  navigation.4  In  this 
system  the  x-axis  points  to  the  polar  North 
(parallel  to  the  LTP)  while  the  z-axis  points  downward  (nadir),  towards  the  Earth’s  surface.  The  y-axis 
completes  the  right-handed  frame  and  points  East  on  the  LTP  (figure  3).  Note  that  while  the  center  of  the 
NED  frame  is  dependent  on  the  body’s  location  (relative  to  ECEF)  the  frame  is  not  a  body- fixed  frame. 

1.  Gravitational  Model 

In  order  to  have  the  body  follow  an  orbital  trajectory,  a  gravitational  model  is  used  to  determine  the  force 
acting  on  the  satellite  at  any  point  in  time.  The  model  implements  the  mathematical  representation  of  the 
geocentric  equipotential  ellipsoid  described  by  the  1984  World  Geodetic  System  (WGS84).2  The  model  is 
able  to  determine  the  Earth’s  gravity  at  any  specific  location  on  the  ECI  frame  limited  only  by  the  geodetic 
height.  By  initializing  latitude,  longitude,  and  altitude  parameters  at  the  expected  launch  date,  any  desired 
orbit  can  be  simulated.  An  advantage  of  utilizing  this  model  in  conjunction  with  6  DOF  dynamics  is  that 


Figure  3:  ECI,  ECEF,  and  NED  reference  frames.  The 
latitude  (A)  and  longitude  (</>)  convention  for  an  arbi¬ 
trary  point  over  the  Earth’s  surface  with  respect  to  the 
ECEF  frame  is  shown. 
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the  satellite  is  not  constrained  to  a  single  or  set  of  predetermined  orbits  but  instead  orbital  transfers  and 
maneuvers  can  be  modeled  together  with  the  attitude. 

2.  Atmospheric  Model 

At  LEO,  the  atmospheric  density  and  composition  have  a  direct  impact  on  the  satellite’s  attitude  especially 
when  one  considers  its  low  mass  and  the  size  of  the  solar  panels  relative  to  the  body.  Using  the  MSISE90 
atmospheric  model  it  is  possible  to  estimate  the  density,  temperature,  and  composition  of  the  atmosphere 
as  well  as  the  number  densities  of  its  components  at  a  specified  altitude  on  the  ECEF.5  Data  from  the 
NRLMSIS-00  (an  updated  version  of  MSISE90)  is  also  used  since  it  contains  additional  data  on  Oxygen 
particles  at  altitudes  above  500  km.6 

3.  World  Magnetic  Model 

The  last  component  of  the  environmental  model  is  a  representation  of  the  Earth’s  magnetic  field  given  by 
the  2010-2015  World  Magnetic  Model  (WMM).  This  model  provides  magnetic  intensity,  inclination,  and 
declination  information  and  a  complete  geometry  of  the  field  in  a  -1  km  to  850  km  range  in  the  NED  frame.7 

C.  External  Disturbance  Models 

For  a  typical  satellite  in  LEO,  there  are  four  significant  attitude  disturbances:  aerodynamic  drag,  magnetic 
field  toques,  gravity-gradient  torques,  and  impingement  by  solar  radiation.8  These  are  regarded  as  separate 
pure  moments  which  correspond  to  the  disturbance  term  in  Eq.  (2). 

1.  Aerodynamic  Drag 

The  aerodynamic  torques  are  produced  by  the  atmospheric  particles  colliding  with  the  satellite  surface. 
Collisions  occur  at  a  higher  frequency  during  maximum  solar  activity  resulting  in  larger  disturbances.  This 
worst  case  scenario  has  been  assumed  for  all  aerodynamic  calculations.  The  aerodynamic  torques  on  the 
satellite  are  estimated  by  Eq.  (10) 


Taero  =  ^pV^CmSl  (10) 

At  a  nominal  altitude  of  500  km  the  NRLMSISE00  model  gives  an  atmospheric  composition  of  94%  O  and 
6%  N.  The  number  and  mass  densities  are  n  =  3.769  x  1014  ks/m-3  and  p  =  1.02  x  1011  ks/m-3  respectively 
and  the  temperature  is  1491K.  The  mean  free  path  between  particles  A  =  27.33  km  which  compared  with  a 
reference  length  l  of  0.1  m  results  in  a  Knudsen  number,  Kn  =  A//  =  273,300,  indicating  that  the  satellite 
lies  in  the  free  molecular  flow  regime. 

Given  the  previous  statement,  the  moment  coefficients  of  the  satellite  were  calculated  using  the  direct 
solution  Monte  Carlo  (DSMC)  program  DAC97,  which  employs  algorithms  based  on  the  methods  described 
by  G. Bird. 9  A  variable  hard  sphere  model  (VHS)  has  been  assumed  for  the  collisions  between  the  particles 
(O  and  N)  and  the  satellite.  V  has  been  assumed  to  be  7.612  km/s,  which  is  equivalent  to  the  mean  orbital 
speed  at  a  500  km  circular  orbit. 

A  total  of  1369  runs  of  the  DAC97  code  have  been  performed  on  the  ARAPAIMA  body  geometry  for  37 
different  a  and  (3  values.  The  angle  of  attack  was  varied  between  -90°  and  90°  in  steps  of  5°;  the  sideslip 
angle  was  similarly  varied  between  0°  and  180°  (the  angle  of  attack  and  the  sideslip  angle  are  defined  as 
seen  in  figure  4).  The  DSMC  results  are  used  to  form  2D  look-up  tables  which  map  the  drag  and  moment 
coefficients  as  functions  of  a  and  /?  as  shown  in  figure  5.  Together  with  the  total  mass  density  obtained  from 
the  MSISE90  model  we  can  use  Eq.  (10)  to  calculate  the  aerodynamic  disturbance  torques  at  any  orientation 
(figure  6,  subfigure  (a)).  The  maximum  force  coefficient  corresponding  to  the  drag  coefficient  in  general 
aerodynamic  terms  is  approximately  2.0  which  is  within  the  expected  range  (2. 0-2. 2)  for  small  satellites  in 
LEO.8 

2.  Magnetic  Residual  Moment 

The  magnetic  disturbance  torque  has  two  major  sources:  the  force  produced  on  a  point  charge  by  the 
magnetic  component  of  the  Lorentz  force,  and  the  torque  experienced  by  an  aspheric  paramagnetic  body 
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Figure  4:  Angle  of  attack  (a)  and  sideslip  angle  (/ 3 )  definitions. 


Aerodynamic  Coefficients 


Figure  5:  Force  and  moment  coefficient  of  the  ARAPAIMA  satellite  as  functions  of  its  orientation  in  terms 
of  a  and  /3. 
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which,  in  the  absence  of  other  torques,  aligns  its  long  axis  with  the  local  magnetic  field.  This  means  that 
all  electrically  conducting  parts  of  the  satellite  contribute  charges  and  produce  time  varying  magnetic  fields. 
Due  to  the  complex  nature  of  evaluating  these  disturbance  torques,  we  resort  to  empirical  data  in  order 
to  estimate  their  effect.  An  estimate  of  the  maximum  magnetic  torque  can  be  obtained  by  combining  all 
contributing  magnetic  effects  into  a  residual  dipole  moment  specific  to  the  satellite  body  and  exposing  it  to 
the  environmental  magnetic  field.8, 10 


Trrnm  —  TTlrrmn  X  B  jt;  (H) 

For  the  ARAPAIMA  satellite,  mrrrirri  was  approximated  to  be  0.1  Am2  aligned  with  the  body  y-axis 
thus  mrrnrn  =  [0,  0.1,0]  Am2.  Together  with  the  WMM  Mrrnrri  can  be  calculated  as  shown  in  figure  6, 
subfigure  (b). 

3.  Gravity  Gradient  Torque 

When  a  satellite’s  center  of  mass  does  not  coincide  with  its  center  of  gravity,  the  variation  of  the  Earth’s 
gravitational  field  over  the  volume  of  the  spacecraft  produces  torques  which  in  the  absence  of  other  distur¬ 
bances  will  try  to  align  one  of  the  body’s  principal  axes  with  the  local  gravity  field  vector.11  This  disturbance 
is  easily  calculated  by 

M99  =  x  Jr  (12) 

The  gravity  gradient  disturbance  contribution  to  the  attitude  dynamics  can  be  calculated  by  evaluating 
Eq.  (12)  along  the  satellite’s  trajectory.  The  resulting  moments  are  shown  in  figure  6,  subfigure  (c). 


Figure  6:  External  disturbance  torques  acting  on  satellite  body  for  a  1  orbit  period 
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4-  Solar  Radiation  Pressure 


When  solar  photons  collide  with  a  satellite,  the  momentum  exchange  between  sunlight  and  the  body  surface 
results  in  a  net  pressure  force  called  solar  radiation  pressure.  This  pressure  is  a  complex  function  of  the 
shape  and  optical  properties  of  the  satellite  as  well  as  the  shading  and  solar  intensity  <f>.  The  worst  case 
scenario  is  given  by  Eq.  (13)  where  is  taken  from  an  STK  simulation  for  the  day  of  January  1st,  2015. 

_*fasa(i  +  z)i, 

IVlgrprp  - 

C 

The  resulting  moment  is  shown  in  figure  6,  subfigure  (d). 

D.  Internal  Disturbance  Models 

External  factors  are  not  the  only  sources  of  disturbances  for  a  satellite;  moving  parts  and  other  mechanical 
interfaces  also  produce  undesirable  torques.  Many  of  these  internal  disturbances  are  also  produced  or  exac¬ 
erbated  as  a  result  of  active  control.  Because  there  is  no  ideal  way  to  mitigate  them,  the  controller  must  be 
able  must  be  able  to  minimize  both  external  and  internal  torques. 


(13) 


1.  Solar  Panel  Flexible  Modes 

In  order  to  deploy  properly,  the  solar  panel  hinges  are  loaded  with  springs  which  rotate  upon  release.  In 
order  for  the  rigid  body  dynamics  outlined  in  section  A  to  be  valid  the  flexible  joints  must  be  modeled 
such  that  any  moment  or  force  resulting  from  their  movement  can  be  folded  into  the  pre-existing  dynamic 
equations.  If  we  regard  the  solar  panels  as  rigid  bodies  the  joints  can  be  modeled  as  a  torsion  spring  by 

Mp  =  JpOp  +  bOp  +  kOp  (14) 

The  resulting  moment  Mp  can  now  simply  be  added  to  Eq.  (2)  as  a  disturbance.  Although  small  in 
comparison  to  the  external  disturbances  (figure  7),  modeling  this  reaction  is  important  in  order  to  monitor 
and  avoid  exciting  structural  modes  through  control  actuation. 


in-ie  Flex  Solar  Panel  Moment 


Time  (s) 

Figure  7:  Flexible  solar  panel  disturbance  moment. 
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2.  Propellant  Slosh 


Even  in  microgravity,  liquid  propellant  slosh  contributes  directly  to  attitude  destabilization.  In  the  case  of 
the  ARAPAIMA  mission,  slosh  can  hinder  mission  completion  by  introducing  unwanted  angular  momentum 
which  could  result  in  a  deviation  from  trajectory  during  orbital  maneuvers  or  in  oscillations  of  the  pointing 
axis  which  prevent  sufficiently  stable  imaging.  For  this  reason,  a  simple  yet  fairly  accurate  spring-mass- 
damper  equivalent  mechanical  model  was  incorporated  in  the  simulation.  This  model  takes  into  account  the 
propellant’s  behavior  through  the  Weber,  Froude  and  Bond  numbers  and  results  in  the  following  moment 
equation  proposed  by  Dodge12  in  terms  of  angular  (0si)  and  linear  (xsi)  displacements 


M-sl  —  (eCc,0  T  rWlsl,oh()')@sl  ^rrflsl,nh'sl,n(%sl,n  T  Hsl,n@sl, o)  T  lsl,n^sl,n-  (15) 

Note  that,  due  to  the  nature  of  the  mass-spring-damper  representation  Eq.  (15),  needs  to  be  modified 
according  to  the  expected  motion  along  each  axis  (especially  depending  on  the  alignment  of  the  thrusters). 
The  x-axis  implementation  is  given  above.  The  output  on  all  three  axes  for  one  orbit  can  be  seen  in  figure  8. 


Figure  8:  Propellant  slosh  disturbance  moment. 


E.  Actuator  Model 

Originally  the  ARAPAIMA  mission  relied  on  reaction  wheels  to  perform  attitude  control,  however,  these 
were  still  dependent  on  the  RCS  thrusters  in  order  to  periodically  despin  the  flywheels.  This  dependency 
coupled  with  volume,  mass,  and  power  constraints  led  to  the  decision  of  eliminating  the  reaction  wheels  in 
favor  of  an  RCS  thruster  actuated  control  system.  A  trade  study  focusing  on  propellant  consumption  and 
attitude  performance  was  made  to  ensure  the  change  would  not  affect  the  satellite’s  ability  to  complete  the 
mission. 

1.  Reaction  Wheels 

The  reaction  wheels  (RWs)  were  modeled  as  brushed  DC  motors  with  an  open-loop  transfer  function  given 

by 


L0rw  -Am  A 

V  (RmJrw)s  T  R-rn  T S  1 
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Adding  a  simple  PI  controller  yields  the  following  closed  loop  transfer  function. 

O  Kh- 

*Lrw  -L^rvi,rw 

T s  H-  (^K^kprw  H-  l)«s  T  Kki  rw 

Acknowledging  that  the  reaction  wheel  angular  momentum  can  be  calculated  as 


H‘ V UJ  J 1 


rw^rw 


then  the  torque  applied  to  the  satellite  is  given  as  a  function  of  its  angular  velocity 


(17) 


Mru)  =  ft  x  Hrw  (18) 

The  propellant  usage  needed  to  offload  the  reaction  wheels  is  determined  by  calculating  the  required 
thruster  burn  time  and  the  saturation  time  of  each  reaction  wheel.  The  burn  time  is  dependent  on  the  mass 
flow  rate  mprop  and  specific  impulse  Isp  of  the  propellant.  The  saturation  time  is  dependent  on  the  reaction 
wheel  nominal  angular  momentum  and  the  root  square  sum  of  the  disturbance  torques  acting  on  the  satellite 
body. 


tburn  —  r.rpi  (19) 

Zil  l^cs 

T 

TTlrprop  =  y  (20) 

J-sp9e 

Thus  the  amount  of  propellant  used  per  offload  maneuver  is  calculated  by: 

TYl b,rw  —  VTlproptb  urn 

For  a  disturbance  scenario  corresponding  to  that  shown  in  figure  6,  evaluating  Eq.  (21) 
a  total  of  0.945  kg  of  propellant  per  orbit  needed  to  offload  the  reaction  wheels.  Note 
case  scenario  which  assumes  constant  imaging  mode  pointing  requirements  (1  aremin). 


(21) 

over  time  estimates 
that  this  is  a  worse 


2.  RCS  Thrusters 

Using  pulse  width  modulation  techniques 
the  thrusters  can  be  “throttled”  at  1% 
increments  of  the  maximum  thrust  value 
over  a  linear  region  made  possible  by 
rapid  solenoid  valve  actuation.  When  op¬ 
erated  within  a  range  of  0-60  Hz,  the 
mass  flow  rate  remains  linear  through  the 
entire  duty  interval  (figure  9).  When  op¬ 
erated  at  higher  frequencies  the  linear  re¬ 
gion  begins  to  deteriorate;  the  current 
model  assumes  a  frcs  =  10  Hz  operating 
frequency. 

Building  on  these  properties  the  RCS 
thruster  model  assumes  a  quantized  lin¬ 
ear  relationship  between  the  commanded 
and  produced  thrust  where  each  step  size 
corresponds  to  a  1%  of  the  maximum 
thrust  value.  A  graphical  representation 
of  this  assumption  is  shown  in  figure  10. 

This  quantized  model  also  introduces 
a  delay  which  corresponds  to  the  differ¬ 
ence  in  sampling  time  of  the  controller, 
valves,  and  onboard  computer.  Because 
the  current  simulation  assumes  the  valves 


Thruster  Valve  Performance  (Frequency  vs  Flow  Rate) 
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Figure  9:  RCS  valves  mass  flow  rate  performance  at  different 
operation  frequencies. 


operate  at  a  slower  frequency  than  the  controller,  this  delay  is 
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currently  equivalent  to  .  Equation  (22)  describes  the  relationship  between  the  commanded  and  delayed 
actuated  moment. 


-M-cmd,(t)  2/rcs^cmd,(t) 

-M-act,(t-\- 1)  ^rcs'-^'act,(t- (-1)  2lrcsTCmd,(t) 


(22) 


Propellant  Consumption  (RCS) 


Figure  10:  Quantization  and  saturation  of  thrust  Figure  11:  Propellant  usage  for  RCS  thruster  dis- 
command  values.  turbance  rejection. 

As  a  means  of  direct  performance  comparison  with  the  reaction  wheel  actuated  control,  the  same  distur¬ 
bance  scenario  was  presented  to  the  thruster  actuated  system.  The  total  propellant  per  orbit  needed  to  reject 
all  exterior  disturbances  was  approximately  0.95  kg  as  shown  in  figure  11.  Figure  12  compares  the  pointing 
performance  of  both  actuators.  Given  that  the  pointing  performance  of  both  actuators  is  comparable  and 
the  propellant  usage  is  the  same  for  this  combination  of  thrust  and  propellant  specifications,  the  RWs  were 
removed  from  the  design  in  order  to  save  mass  and  volume  needed  for  other  components. 


Figure  12:  RW  and  RCS  thruster  point  performance 
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III.  Controller  Design 


The  ARAPAIMA  mission  requires  that  its  ADCS  be  able  to  handle  long  angle  transition  maneuvers 
as  well  as  precise  attitude  tracking  imaging  maneuvers.  The  approach  presented  here  uses  a  combination 
of  nonlinear  eigenaxis  control  and  PID  controllers  about  each  axis.  For  the  moment,  a  periodic  switching 
scheme  is  presented  together  with  a  bump  conditioning  algorithm  with  scheduled  gains  which  attempts  to 
smooth  the  transition  between  controllers. 

A.  Eigenaxis  Control 

The  large  angle  maneuver  control  law  is  based  on  the  quaternion  feedback  regulator  proposed  by  Wie13  and 
applied  by  Pong14  to  CubeSat  missions.  The  controller  consists  of  a  linear  error  quaternion  feedback  and  both 
linear  and  nonlinear  angular  rate  feedback  which  counteract  gyroscopic  coupling  torques.13  This  approach 
is  based  on  eigenaxis  rotations  and  is  analogous  to  the  well  known  PD  controller.  It  was  designed  with  large 
angle  maneuvers  in  mind  and  its  global  stability  as  well  as  its  robustness  to  inertia  matrix  uncertainty  has 
been  proven  in  relation  to  various  spacecraft  applications.13, 15,16  The  torque  command  given  by  the  control 
law  is  expressed  in  the  Laplace  domain  as 


Mcrnd  =  -a?  x  Ju  -  Kdu;  -  Kpqe  (23) 

where  the  proportional  and  derivative  gain  matrices  are  given  by  Kp  =  pJ  and  Kd  =  dJ  respectively;  p 
and  d  are  scalar  tuning  parameters. 

B.  PID  Control 

As  with  all  PD  controllers,  the  lack  of  an  integral  term  introduces  non-zero  steady  state  errors.  In  order 
to  achieve  the  high  accuracy  needed  for  imaging  modes,  a  simple  PID  controller  is  applied  about  each  axis. 
Taking  the  vector  part  of  the  error  quaternion 


MCmd  —  —  [Kp  -f — - — j-  Kds\qev  (24) 

where  the  gains  again  are  the  product  of  the  body  moment  of  inertia  matrix  and  scalar  tuning  parameters. 

However,  high  precision  integral  controllers  are  ill  suited  for  large  angle  maneuvers  and  are  heavily 
dependent  on  the  magnitude  of  the  error  or  initial  conditions  to  ensure  stability.  Therefore,  the  PID  control 
law  is  implemented  as  an  extension  of  the  eigenaxis  controller  which  comes  into  effect  once  a  certain  set  of 
circumstances  are  met.  These  circumstances  are  a  weighed  combination  of  the  current  attitude  error,  its 
behavior  over  time,  and  the  current  mission  operational  mode.  This  way,  the  eigenaxis  controller  operates 
throughout  the  majority  of  the  maneuvers  and  ensures  the  appropriate  conditions  are  met  for  the  PID 
controller  to  take  over.  If  the  error  grows  past  a  predetermined  set  of  bounds,  operation  will  switch  back  to 
the  eigenaxis  controller  to  prevent  instability. 

1.  Bump  Conditioning 

The  smoothness  with  which  the  switch  or  transfer  between  controllers  occurs  is  addressed  by  considering 
bump  conditioning  techniques  as  presented  by  Peng.17  Although  complex  system  dynamics  prevent  com¬ 
pletely  bumpless  transitions,  some  improvement  is  achieved  by  altering  Eq.  (24)  to  the  following  form: 

MCmd,(t-\-l)  —  —KpQev,(t)  —  (Kt{MCmd,(t)  ^ cmd,(t ))  +  Qev,(t))  ~  ^dsQev,(t)  (25) 

Note  that  the  difference  between  Eq.  (24)  and  Eq.  (25)  lies  on  how  the  integral  action  is  addressed. 
Because  the  eigenaxis  controller  lacks  integral  terms,  this  kind  of  conditioning  in  only  useful  when  switching 
to  PID  control;  similar  corrections  are  unnecessary  when  operation  switches  back  to  eigenaxis  control. 

The  effectiveness  of  Eq.  (25)  can  be  augmented  by  the  use  of  scheduled  gains.  Building  on  the  notion 
that  integral  controllers  are  sensitive  to  abrupt  changes,  linear  transitions  can  be  scheduled  such  that  the 
gain  values  are  identical  to  those  of  the  eigenaxis  controller  at  the  instant  when  the  switch  to  PID  occurs. 
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This  further  smooths  the  introduction  of  integral  action  into  the  system  since  Ki  would  start  at  zero  and 
then  increase  according  to  some  liner  function  of  time  f(t)  as: 

P2  =Pl+f(t)(p2-Pl) 

h  =  h+  -  h)  (26) 

d2  ~  di  + f(t)(d2  -  dx) 

where  /(£)  ranges  from  0  to  1.  The  slope  of  f(t)  then  controls  the  abruptness  of  the  transition. 

IV.  Numerical  Simulation  and  Results 

A  numerical  simulation  model  has  been  prepared  in  a  combined  MATLAB  and  Simulink  environment. 
This  setup  has  the  advantage  of  including  many  blocksets  and  functions  which  simplify  the  calculations 
needed,  particularly  the  built-in  interface  for  the  environmental  models  mentioned  in  section  B.  The  model 
runtime  was  set  to  simulate  a  single  92  min  orbit  in  close  proximity  (250  m)  to  the  RSO  using  a  4th  order 
fixed-step  extrapolation  integrator  (odel4x)  with  a  step  size  of  0.05  s.  The  RSO  and  satellite  orbits  were 
determined  from  STK  simulations  used  to  design  the  ARAPAIMA  mission  trajectories.  These  give  us  the 
required  c ccrrid  and  qCTnd  commands  needed  to  track  the  RSO  (i.e.  controller  commands). 

The  model  takes  into  account  all  dynamics  outlined  in  section  II  and  applies  both  control  laws  discussed 
in  section  III.  However,  due  to  the  short  span  of  the  simulation,  the  switching  between  controllers  has  been 
made  periodic  in  order  to  show  the  characteristics  of  both  control  laws.  Figure  13  shows  the  combined 
internal  and  external  disturbance  torques  acting  on  the  body  throughout  the  simulation.  For  the  purposes  of 
initializing  all  environmental  models,  the  date  is  assumed  to  be  January  1st  2015.  This  date  also  corresponds 
to  the  date  of  the  RSO  orbit  simulation. 


Figure  13:  Flexible  solar  panel  disturbance  moment. 


The  eigenaxis  controller  was  tuned  so  that  the  scalar  tuning  parameters  pi  =  8,  d\  =0,  and  i\  =  0  whereas 
the  PID  controller  was  tuned  so  that  P2  =  2,  d 2  =  15,  and  i 2  =  0.1. Figure  14  shows  the  overall  pointing 
performance  of  the  satellite  for  different  control  settings.  In  order  to  assess  the  advantages  of  switching 
while  still  quantifying  the  smoothness  of  the  transition  between  controllers  three  different  scenarios  are 
examined.  The  first  (figure  14,  sub  figure  (a))  shows  the  pointing  performance  with  no  bump  conditioning 
or  gain  scheduling.  This  gives  a  steady  state  error  of  0.18  arcmin  for  the  eigenaxis  controller  and  0.025  for 
the  PID  controller.  However,  the  switch  causes  significant  overshoot  (1-7  arcmin).  Note  that  the  PID  is 
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Time  (s) 

(a)  Pointing  performance  with  no  bump 
conditioning  and  no  gain  scheduling 


Time  (s) 

(c)  Pointing  performance  with  bump  con¬ 
ditioning  and  gain  scheduling 


Pointing  Error  With  Bump  Conditioning 


(b)  Pointing  performance  with  bump  con¬ 
ditioning  but  no  gain  scheduling 


Switch  Performance  Comparison 


(d)  Bump  conditioning  performance  com¬ 
parison 


Figure  14:  Switching  controller  performance. 


sensitive  to  changes  in  angular  acceleration  whereas  the  eigenaxis  controller  remains  at  steady  state  except 
for  short  jumps  where  the  accuracy  improves. 

The  second  scenario  (sub figure  (b))  shows  the  performance  when  the  PID  controller  is  as  in  Eq.  (25). 
This  offers  noticeable  improvement  in  the  switch  overshoot  (0.65-3  arcmin)  but  because  the  Kt  term  remains 
non-zero  after  the  transient  period,  it  offers  no  improvement  on  the  overall  performance.  This  is  fixed  in  the 
third  scenario  when  gain  scheduling  is  introduced.  The  scheduled  gains  assume  a  transient  period  of  200  s 
during  which  the  Kp ,  Ki ,  Kd  gains  are  changed  according  to  Eq.  (26)  and  are  shown  in  figure  15.  The  bump 
conditioning  gain  does  not  behave  this  way,  but  instead  is  a  fixed  non-zero  value  during  the  transient  period 
and  is  zero  elsewhere.  This  approach  further  improves  the  overshoot  (0.70-1.20  arcmin)  and  restores  the 
improvement  on  pointing  performance  seen  in  the  first  scenario  yet  at  the  cost  of  longer  transition  periods. 
Subfigure  (d)  gives  a  side  by  side  comparison  of  the  bump  conditioning  given  by  each  approach. 
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Gain  Scheduling  Transients 


Figure  15:  Gain  scheduling  transients  corresponding  to  those  in 
figure  14,  subfigure  (d) 


A.  Future  Work 

While  the  current  results  are  promising,  there  remain  ways  in  which  the  fidelity  of  the  model  might  be 
improved  in  order  to  further  test  and  prove  the  stability  of  the  ADCS.  As  the  control  laws  must  be  imple¬ 
mented  on  the  satellite  on-board  computer,  the  entire  model  must  be  discretized  and  made  to  account  for 
the  different  sampling  rates  and  inherent  delays  of  digital  implementation.  In  direct  relation  is  the  matter  of 
hardware  noise  and  interference  along  with  a  corresponding  filter  design  to  compensate  for  the  rate  transition 
and  inaccuracies  of  sensor  readings.  Figure  16  gives  an  example  of  controller  performance  with  unfiltered 
sensor  noise;  the  noise  levels  are  based  on  the  real  mission  hardware. 


Pointing  Performance  With  Noise  Pointing  Performance  With  Noise  (4500  s  Window) 


Time  (s)  Time  (s) 

Figure  16:  Pointing  error  with  sensor  noise  in  the  system. 


It  is  clearly  shown  how  the  sensor  noise  has  a  significant  impact  on  pointing  accuracy  especially  while 
during  eigenaxis  control  since  from  Eq.  (23)  we  can  see  that  it  is  directly  susceptible  to  both  STR  and  IMU 
noise.  The  PID  maintains  an  accuracy  of  0.5  arcmin  but  in  a  real  scenario  an  IMU  samples  at  a  much  faster 
rate  than  a  STR  so  attitude  estimate  propagation  would  still  be  affected  by  IMU  noise,  decreasing  the  PID 
performance. 
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V.  Conclusion 


The  purpose  of  this  paper  has  been  to  demonstrate  the  validity  of  the  ADCS  design  approach  as  well  as 
the  capability  of  CubeSat  class  satellites  to  contribute  to  current  areas  of  interest  using  COTS  products  and 
without  the  necessity  of  a  multi- million  dollar  budget.  Going  further  than  a  feasibility  analysis  the  aim  is 
to  prove  through  testing  and  design  the  performance  of  a  real  system  that  is  able  to  extend  into  hardware 
integration  and  implementation.  Through  the  use  of  proven  models  and  methods  and  a  combination  thereof 
it  has  been  shown  that  the  hardware  (through  specifications)  and  design  are  able  to  achieve  the  performance 
stipulated  in  the  mission  criteria.  Future  work  and  necessary  improvements  are  expected  to  increase  the 
quality  and  reliability  of  the  working  system. 
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